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Abstract 

A  method  is  proposed  for  super-resolving  multi-channel  data  with  applications  to 
PREDATOR  video  sequences.  Using  a  generalization  of  Papoulis’  sampling  theorem, 
a  closed-form  solution  has  been  obtained  leading  to  a  high-speed  algorithm  which  can 
be  realistically  applied  to  large  data  sets  such  as  video  sequences.  In  existing  multi¬ 
frame  methods  it  is  a  common  practice  to  assume  that  the  channel  transfer  functions 
are  known  and  invariant  from  one  frame  to  another,  using  empirical  models  such  as 
Gaussian,  sine,  etc.  We  have  assumed  that  the  transfer  functions  are  unknown  and 
may  vary  even  when  the  same  sensor  is  employed,  and  hence  use  the  observed  data 
to  derive  the  Point  Spread  Function  (PSF)  for  each  frame.  The  estimated  PSFs  are 
used  in  the  super-resolution  algorithm.  Results  on  PREDATOR  video  images  are 
then  given. 

Keywords:  Super-resolution,  video  data  processing,  multi-channel  sampling,  data- 
driven  estimation 
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1  Introduction 


Super-resolution  refers  to  methods  that  attempt  to  increase  resolving  power  by  means 
of  image  processing.  Super-resolving  video  data  in  particular  is  potentially  of  signifi¬ 
cant  interest  in  visual  surveillance  and  monitoring  of  human  and  vehicular  activities 
in  both  civilian  and  battlefield  applications  [10]  [23].  In  fact,  many  automatic  target 
recognition,  detection  and  identification  problems  suffer  from  lack  of  adequate  reso¬ 
lution  of  the  data  (PREDATOR  video  is  a  good  example).  When  replacing  a  sensor 
is  not  a  feasible  solution  (due  to  either  cost  or  technical  limitations),  post-processing 
by  means  of  super-resolution  algorithms  is  the  only  means  of  improving  the  resolving 
power  for  enhancing  performance. 

Earlier  work  on  super-resolution  involved  single-frame  methods  in  which,  assuming 
that  the  signal  was  of  compact  support  in  space  (time),  its  spectrum  was  extrapo¬ 
lated  outside  a  known  support.  Variations  of  single-frame  methods  include:  methods 
based  on  analytic  continuation  [5],  those  based  on  prolate  spheroidal  wave  functions 
[20][21][22],  the  well-known  Papoulis-Gerchberg  [3][12]  algorithm,  and  some  proba¬ 
bilistic  methods  [7]. 

Recent  work  has  mostly  involved  multi-frame  data  [1]  [6]  [8]  [9]  [17]  [19],  where  the 
assumption  is  that  the  signal  is  not  space-limited  but  rather  band-limited  and  that 
each  frame  has  been  down-sampled  so  that  every  frame  is  an  aliased  representation 
of  the  underlying  signal.  By  making  the  assumption  that  the  sampling  matrices  are 
not  correlated  one  can  then  attempt  to  increase  the  sampling  density  by  trading  off 
the  temporal  bandwidth,  and  hence  to  unscramble  some  or  all  of  the  aliased  portion 
of  the  spectrum. 

In  this  paper  we  will  be  mainly  concerned  with  multi-frame  methods,  which  clearly 
are  better  adapted  to  processing  video  data.  However,  existing  methods  suffer  from 
a  number  of  shortcomings  which  we  have  addressed: 

•  In  the  existing  literature,  the  sensor  PSF  is  assumed  invariant  from  one  frame 
to  another;  empirical  models  such  as  Gaussian,  sine,  etc.  are  generally  used. 
The  analysis  is  thus  simplified  at  the  cost  of  reduced  performance.  We  propose 
a  method  where  the  PSF  is  adaptively  estimated  from  the  observed  data  for 
each  frame  in  the  sequence.  By  allowing  different  PSFs  for  different  frames,  our 
method  extends  multi-frame  methods  to  a  true  multi-channel  framework. 

•  Speed  is  an  important  issue  when  dealing  with  video  sequences.  Most  existing 
methods  are  slow.  We  have  proposed  a  closed-form  solution,  resulting  in  a 
high-speed  algorithm  and  allowing  for  realistic  application  to  video  data. 
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2  Multi-channel  super-resolution 


Our  problem  is  formulated  as  follows:  We  have  a  sequence  of  M  images  that  were 
collected  by  a  set  of  systems  whose  transfer  functions  are  unknown.  The  output  of 
each  channel  is  down-sampled,  in  the  sense  that  each  image  alone  cannot  represent 
the  underlying  input  signal.  We  are  interested  in  establishing  a  method  of  recon¬ 
structing  the  input  signal,  as  well  as  conditions  under  which  the  reconstruction  leads 
to  a  physically  meaningful  solution. 

Our  analysis  is  based  on  a  generalization  of  Papoulis’  sampling  theorem  [13].  Con¬ 
sider  the  very  simple  example  of  reconstructing  a  ID  signal  observed  by  two  identical 
sensors  sampling  uniformly  with  the  same  period  T  but  with  a  temporal  shift  of  St 
between  them.  If  St  is  an  integer  multiple  of  T  then  the  two  sensors  will  produce 
redundant  information,  using  twice  as  much  temporal  bandwidth  as  required.  There¬ 
fore,  in  order  to  be  able  to  trade  off  some  of  the  temporal  bandwidth,  we  must  have 
St  ^  kT  'ik  G  Z.  In  general,  this  sampling  is  non-uniform  unless  St  =  {k  +  |)r 
for  some  integer  k.  However,  the  samples  constitute  a  special  case  of  non-uniform 
sampling  referred  to  as  non-uniform  recurring  samples,  bunched  samples  or  interlaced 
samples.  This  is  shown  in  Figure  1. 


Figure  1:  Multi-channel  interlaced  sampling 

Therefore,  even  in  the  simplest  case,  the  problem  is  outside  the  scope  of  uniform 
sampling  theory.  In  fact,  in  practice  the  problem  can  be  much  more  complicated  due 
to  the  following  factors: 

•  In  higher  dimensions  inter-frame  shifts  will  generalize  to  other  transformations 
whose  parameters  are  in  general  unknown  and  need  to  be  accurately  estimated 
(the  images  need  to  be  registered  at  sub-pixel  accuracy). 

•  Sensors  in  general  cause  degradations  due  to  their  PSFs. 
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•  Sensors  may  be  different  (have  different  PSFs),  or  in  the  case  when  the  same 
sensor  is  used,  the  PSF  for  each  may  vary  frame  due  to  lighting  conditions, 
noise,  etc. 

Below,  we  will  first  set  the  theoretical  background  for  designing  our  algorithm. 


2.1  Theory 


Definition  2.1  A  sttll{C)  C  H”  is  said  to  be  a  set  of  uniqueness  for  a  given  class 
of  functions  C  if  f  =  0  on  U{C)  implies  that  f  is  the  null  function:  /  =  0. 

In  general,  a  set  of  uniqueness  does  not  have  to  be  discrete,  but  all  sets  of  uniqueness 
considered  here  will  be  discrete  subsets  of  H",  i.e.  multivariate  sequences.  Also 
the  functions  we  will  be  dealing  with  will  belong  to  the  class  of  entire  functions  of 
exponential  type,  namely  functions  that  belong  to  a  Paley- Wiener  space  V  [11]  over 
a  compact  support  B  C  IR”.  The  following  theorem  will  then  generalize  Papoulis’ 
sampling  theorem  to  a  multi-dimensional  multi-resolution  framework,  where  multi¬ 
resolution  implies  that,  contrary  to  Papoulis’  theorem,  the  channel  sampling  densities 
do  not  have  to  be  equal,  provided  that  the  overall  density  equals  that  of  Nyquist. 

Theorem  2.1  (Multi-Resolution,  Multi-channel  Band-limited  Sampling)  Let  /(t)  € 
'P(B)  be  convolved  with  a  set  of  mutually  uncorrelated  linear  shift-invariant  filters 
in  whose  spectra  {ftm(u)}m=i,...,M  are  uniformly  bounded  away  from  zero 

in  B.  Assume  also  that  the  outputs  are  known  only  on  a  discrete  set  of  points 
{5'm(Smk:)}m=i,...,Af,  where  the  nonsingular  matrices  S,n  are  such  that 

M 

U{S»l‘}k52.=M{P(B))  (1) 

m=l 

Moreover,  let  {</?m(t)}m=i,...,M  be  a  set  of  linear  shift-invariant  functions  in  4(IR”) 
satisfying  the  following  biorthogonality  conditions  overB; 

=  1b  (2) 

Then  there  exist  scalars  c^k  such  that 

1  M 

/(f)  =  "77  ^  —  Smk)  (3) 

m=l  fceZ" 

with  convergence  being  uniform  and  —  fl'm(Smk). 
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For  a  proof  see  [16]. 

In  fact,  {v?m(t  -  which  we  shall  refer  to  as  the  biorthogonal  set, 

constitutes  a  Riesz  basis  in  V(B),  and  each  component  in  the  set  can  be  found  as  a 
space-varying  linear  transformation  of  Papoulis’  periodic  inverse  filters  (see  [16]  for 
details). 

In  contrast  to  Papoulis’  theorem,  where  the  reconstruction  algorithm  is  based  on 
solving  a  large  linear  matrix  equation,  the  above  theorem  provides  a  more  straight¬ 
forward  method.  In  fact,  note  that  the  biorthogonal  set  is  known  if  the  channel 
transfer  functions  (i.e.,  the  PSFs)  can  be  estimated  and  the  sampling  matrices  can  be 
specified.  In  image  processing  terminology,  the  latter  implies  that  the  image  frames 
need  to  be  registered.  For  the  registration  we  have  used  a  method  based  on  phase 
correlation  developed  in  [2]  [18]. 

As  for  the  reconstruction,  in  principle,  any  basis  can  be  used  (wavelet,  orthogonal 
polynomials,  etc.).  The  Fourier  basis,  for  instance,  can  be  found  using  the  duality 
due  to  the  Paley- Wiener  theorem  [11]: 

{*^m(u)exp(— ?(STOk,  u))}to=i,...,m  (4) 

where,  for  our  application,  the  k’s  take  their  values  in  compact  subsets  of  {)  is 
the  usual  inner  product  in  IR^,  and  u  (with  abuse  of  notation)  will  hereafter  stand  for 
discrete  frequency  vectors  sampled  at  super-resolution  density.  Note  from  equation 
(4)  that  the  expansion  of  each  frame  on  the  biorthogonal  basis  can  be  found  by 
expanding  on  the  standard  orthonormal  basis  followed  by  applying  the  biorthogonal 
filter  :,3to(u)  (a  projection).  The  basis  is  completely  specified  if  we  can  find  which 
are  in  turn  specified  according  to  the  biorthogonality  condition  in  terms  of  the  channel 
transfer  functions.  Once  the  channel  transfer  functions  have  been  estimated,  we  can 
reconstruct  the  biorthogonal  set  (or  a  pseudo-biorthogonal  set  if  the  transfer  functions 
contain  singularities).  Therefore,  in  the  next  section  we  will  develop  a  data-driven 
method  for  estimating  the  PSFs. 

2.2  Estimating  the  PSFs 


In  this  section  we  will  define  our  PSF  model,  or  more  precisely  its  spectrum,  ref- 
ered  to  as  the  Optical  Transfer  Function  (OTF).  We  will  then  develop  a  method 
for  data-driven  estimation  of  the  OTFs.  As  we  will  see  in  the  next  section,  the 
proposed  method  is  particularly  applicable  to  video  data,  since  it  requires  two  im¬ 
ages  with  only  relative  shifts  between  them.  In  fact,  in  video  sequences  (such  as 
PREDATOR  data)  transformations  between  successive  frames  can  be  closely  approx¬ 
imated  by  shifts  and  a  small-angle  rotation  within  the  image  plane.  Therefore,  by 
using  a  registration  method  (see  the  section  on  implementation),  successive  frames 
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are  rotation-compensated  prior  to  PSF  estimation. 

In  the  existing  literature  on  super-resolution,  it  is  a  common  practice  to  assume 
that  the  OTF  is  known  and  traditionally  empirical  models  such  as  Gaussian,  sine, 
etc.  are  adopted.  Although  these  models  can  greatly  reduce  the  complexity  in  both 
analysis  and  implementation  due  to  their  nice  global  behavior  (isotropy  around  the 
origin),  they  fail  to  capture  local  variations  of  the  OTF  and  hence  are  not  useful  for 
adaptive  estimation  and  inversion. 

We  will  show  below  how  more  versatile  local  adaptive  models  can  be  built  using 
local  isotropy.  For  this  purpose  we  will  assume  that  the  OTFs  of  the  imaging  systems 
are  linear  shift-invariant  and  of  finite  duration  and  hence  can  be  represented  by  FIR 
filters: 

^m(u)  =  Cfcexp(-i(u,k))  (5) 

keT 

where  the  c^s  are  constants,  and  T  is  a  compact  subset  of  Z^. 

As  is  well  known,  due  to  the  fundamental  theorem  of  algebra,  in  the  univariate 
case  the  transfer  function  hm  can  be  factored  as  the  product  of  its  roots: 

N 

hm{u)  =  A  JJ(1  -  «fcexp(-iu)) 
fc=0 
N 

=  A'^{1  -  akZ~'^) 

A:=0 

where  N  specifies  the  extent  of  the  impulse  response,  the  a^s  and  A  are  constants, 
and  2  =  exp(iu). 

Clearly  hm  is  then  specified  up  to  a  scale  factor  if  its  roots  are  known.  Unfortu¬ 
nately,  due  to  the  general  impossibility  of  factoring  polynomials  in  higher  dimensions, 
this  simple  convenient  factorization  of  the  spectrum  in  univariate  problems  does  not 
extend  to  n  >  1  dimensions. 

However,  we  will  show  that  in  the  two-dimensional  case  the  problem  can  have  a 
solution  too.  We  start  by  assuming  that  the  zeros  of  our  channel  transfer  functions 
occur  at  isolated  points^ .  Since  the  transfer  functions  are  then  either  locally  convex 
or  concave  around  their  zeros,  one  may  assume  local  isotropy  in  small  neighborhoods 
of  the  roots.  The  following  results  will  show  how  a  two-dimensional  spectrum  with 
isolated  zeros  can  be  factored. 


(6) 

(7) 


^In  general,  the  zeros  of  an  entire  function  of  two  or  more  variables  do  not  have  to  be  isolated. 
However,  in  practice  they  often  occur  at  isolated  points  [15]. 


5 


Proposition  2.1  Let  {Cn}n^z  ®  set  of  isolated  points  in  C".  Then  any  function 
defined  over  C  ”,  which  is  locally  isotropic  around  these  points,  can  be  decomposed 
into  a  sum  of  functions,  each  globally  isotropic  around  one  point  of  the  set. 

The  existence  of  such  a  decomposition  under  appropriate  boundary  conditions  is 
obvious.  However,  the  decomposition  is  not  necessarily  unique,  which  suggests  that 
the  inverse  may  not  necessarily  be  true.  In  other  words,  a  sum  of  functions  that  are 
globally  isotropic  around  distinct  isolated  points  will  not  necessarily  yield  a  function 
that  is  locally  isotropic  around  these  points.  In  fact,  counter-examples  can  be  readily 
found. 


Theorem  2.2  Let  hk{u,  v)  be  the  spectrum  of  a  hand-limited  function  vanishing  out¬ 
side  a  compact  support  T  C  IR^  Let  also  hk{u,v)  be  isotropic  around  the  point 
{fJ'kjVk).  Then 

Hm  bkhk{u,  v)  =  exp{-iu>tk)  (8) 

where  u  =  {{u-  Uk)^  (v  -  VkY)^  and 
ik  =  (xk  +  yl)^  with  {xk,yk)  €  T. 

For  a  proof  see  the  appendix. 

Let  us  now  assume  that  the  channel  transfer  functions  are  locally  isotropic  around 
their  roots.  Then  according  to  Proposition  2.1  we  can  write 

^  A 

hm{u,v)  =  22^mkiu,v)  (9) 

k 

where  for  every  m,  h,nk  is  globally  isotropic  around  one  root  of  hm-  We  shall  denote 
the  roots  of  hm  by  the  set  where  iif  is  a  compact  index  set. 

Therefore,  an  OTF  which  is  locally  isotropic  around  its  roots  can  be  expanded  as 
follows: 

h{u,v)  =  'Y^akexp{—iutk)  where  ak  =  bf^  (10) 

k 

To  physically  interpret  this  convergence  at  the  limit,  notice  that  the  bkS  define  the 
radii  of  local  isotropy  around  the  roots.  Therefore,  the  more  local  the  isotropy  is, 
the  better  the  OTF  can  be  approximated  by  an  expansion  of  the  form  (10).  This 
assumption  of  local  isotropy  clearly  relaxes  the  severe  constraints  usually  imposed  in 
the  literature  by  global  symmetries. 

2.3  Spectral  Factorization 

In  the  classical  literature  this  term  is  used  when  it  is  required  to  find  a  function  whose 
power  spectrum  is  known.  The  problem  has  a  solution  if  the  power  spectrum  satisfies 
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the  Paley- Wiener  condition  [14].  Below  we  will  use  the  term  in  a  more  general  context 
where  two  functions  are  sought  whose  cross  power  spectrum  is  known.  We  will  show 
that  the  solution  can  be  found  if  the  functions  are  of  compact  support  (e.g.  FIR 
filters).  The  support  constraint  condition  is  clearly  equivalent  to  the  Paley- Wiener 
condition  due  to  the  uncertainty  principle  [14]. 

From  the  expansion  in  (10)  and  the  fundamental  theorem  of  algebra,  it  immedi¬ 
ately  follows  that 


A(u,u)  =  A]][(l  -  Cfcz  ^)  (11) 

k 

where  A  and  the  c^s  are  constants  and  2  =  exp(ia;)  =  exp  ((u  -  ukf  +  (v  -  . 

Clearly,  the  factorization  is  only  valid  if  the  c^s  are  identically  equal  to  unity. 
Therefore,  if  the  roots  of  the  OTF  are  specified,  the  OTF  is  known  up  to  a  scale 
factor.  We  will  deal  with  the  scale  factor  shortly.  Let  us  first  see  how  the  roots  of  the 
OTF  can  be  estimated  from  two  observations  of  the  same  scene. 

Consider  the  case  where  an  object  function  has  been  observed  by  two  systems 
modeled  as 

=  hfi  +  «i  (12) 

32  =  ^2/2  +  ^2  (13) 

where  /2  is  a  shifted  version  of  fi,  hi  and  ^2  are  the  OTFs,  ni  and  n2  are  random 
noise  processes,  and  and  32  are  the  observed  image  spectra. 

When  /i  and  /2  are  uniformly  bounded  away  from  zero  on  their  supports,  we  can 
write 

di  =  hgifi  (14) 

^2  =  ^52/2  (13) 

where  hgi  =  hi  +  hif^  and  hg2  =  ^2  +  ^2/25  with  *  denoting  the  complex  conju¬ 
gate.  We  shall  refer  to  such  transfer  functions  as  the  Generalized  Transfer  Functions 
(GTFs)  of  the  imaging  systems. 

We  can  now  note  that  the  magnitude  of  the  cross  power  spectrum  is  given  by 

I  Mi  1=1  I  (16) 

Therefore,  since  both  hgi  and  hg2  are  assumed  to  be  FIR  filters,  the  poles  of  the 
cross  power  spectrum  will  be  uniquely  determined  by  the  roots  of  hg2-  In  practice  this 
implies  that  the  roots  of  hg2  will  be  present  in  the  form  of  singularities,  appearing 
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as  a  set  of  spikes  scaled  at  different  frequencies.  In  other  words,  the  roots  can  be 
identified  by  simple  inspection  of  the  magnitude  of  the  cross  power  spectrum.  Figure 
2  shows  an  example  of  these  spike  patterns. 


Figure  2:  Typical  spike  patterns  appearing  in  the  magnitude  of  a  cross  power  spectrum 
corresponding  to  isolated  zeros  of  one  of  the  transfer  functions 

To  determine  the  scale  value  A  we  will  assume  that  the  GTFs  preserve  the  mean 
value  of  the  object  function  and  hence 


nfc(l  -  ea:p(-i7yfc)) 

where  7^:  =  {ul  + 

Therefore,  by  identifying  the  singularities  (spikes)  in  the  magnitude  of  the  cross 
power  spectrum,  we  can  completely  specify  the  corresponding  GTF. 

2.4  Constructing  Pseudo-biorthogonals 

In  the  absence  of  noise  the  GTF  is  equal  to  the  OTF  and  when  the  OTF  is  an  ideal 
all-pass  filter  it  corresponds  to  the  noise  process.  We  shall  now  use  a  simple  method 
for  constructing  the  biorthogonaJ  set.  In  fact,  since  in  our  model  the  GTFs  (or  the 
OTFs)  are  assumed  to  have  a  set  of  isolated  zeros,  we  can  only  attempt  to  construct 
a  pseudo-biorthogonal  set.  One  may  consider  several  approaches  for  constructing  a 
stable  pseudo-biorthogonal  set,  for  instance  by  using  an  approximation  around  the 
origin. 

Below,  we  will  use  a  simple  method  which  is  based  on  the  following  identity  derived 
from  the  biorthogonality  condition: 

(1  +  0m)~^  =  hmO-  +  hm)~^ 
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(18) 


Using  a  first-order  expansion  of  the  left-hand  side  around  the  origin,  we  find  that 


1  ~  ~  ^m(l  +  hm)  ^ 

and  hence 

h* 

m  ~  _ 122 - 

Ym  —  ^  ^  ; 

h  h*  4-  h* 

which  has  some  resemblance  to  the  standard  Wiener  filter. 


(19) 

(20) 
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3  Implementation  and  Results 


We  have  implemented  our  method  using  a  Fourier  basis.  The  following  are  the 
steps  in  the  algorithm: 

•  For  a  sequence  of  M  successive  frames  the  PSFs  are  estimated  by  estimating 
the  inter-frame  rotations  [2]  and  then  using  the  method  described  above. 

•  The  pseudo-biorthogonals  are  constructed  for  projection  onto  optimal  bases. 

•  Each  frame  is  decomposed  on  a  standard  orthonormal  Fourier  basis  and  then 
projected  onto  a  more  optimal  Riesz  basis  using  the  corresponding  pseudo- 
biorthogonal  transfer  function. 

•  Using  sub-pixel  inter-frame  displacement  values  [18],  successive  frames  are  re¬ 
aligned  and  added  into  a  single  frame  according  to  (3). 

The  method  has  been  tested  on  numerous  images  from  PREDATOR  video  sequences, 
some  of  which  are  shown  below.  The  results  have  been  compared  to  interpolation  by 
Shannon’s  sampling  expansion  (i.e.  infinite-order  interpolation)  followed  by  image 
sharpening. 

In  each  sequence  presented  below,  image  (a)  shows  one  of  the  low-resolution  frames 
in  a  sequence  of  images;  image  (b)  is  its  interpolated  version;  and  image  (c)  is  the 
super-resolved  version.  Interpolating  using  Shannon’s  sampling  expansion  is,  in  fact, 
equivalent  to  decomposing  the  frame  on  the  standard  orthonormal  basis  and  then 
ideally  low-pass  filtering  for  reconstruction  on  a  denser  basis.  Clearly,  lower-order 
interpolations  (such  as  bilinear)  would  yield  even  lower  quality  images.  The  interpo¬ 
lations  in  (b)  have  two  main  shortcomings  compared  to  the  super-resolved  images  in 

(c): 

•  As  is  well  known,  an  orthonormal  basis  is  only  sub-optimal  in  the  presence  of 
noise  and  other  degradations  and  hence  in  practice  Shannon’s  interpolation  can 
only  provide  sub-optimal  results.  Although,  when  the  signal  to  noise  ratio  is 
high,  interpolation  followed  by  sharpening  may  result  in  good  results,  for  a  very 
low  signal  to  noise  ratio  this  approach  can  considerably  enhance  the  noise  and 
the  artifacts  of  the  contaminated  signal.  Therefore,  an  advantage  of  our  super¬ 
resolution  algorithm  is  its  capability  of  resolving  visual  information  in  noisy 
data  such  as  PREDATOR  video. 

•  When  using  Shannon’s  interpolation,  samples  are  only  specified  by  a  single 
frame  and  hence  the  temporal  bandwidth  is  not  exploited,  as  in  our  super¬ 
resolution  method,  for  increasing  the  sampling  density. 
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Note  that  the  images  used  for  experimentation  are  real  data  for  which  the  ground 
truth  is  not  available  and  therefore  it  is  not  possible  to  investigate  the  quality  of 
the  super-resolved  or  interpolated  images  in  terms  of  signal-to-noise  ratio  or  other 
quantifying  measures.  However,  our  future  plan  is  to  investigate  the  quality  of  the 
results  in  terms  of  the  performance  improvement  in  state-of-the  art  automatic  target 
recognition,  detection  and  identification  methods.  Results  should  be  available  in  a 
comprehensive  report  in  the  near  future. 
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In  the  following  images  we  have  zoomed  in  to  some  regions  in  Figure  5  and  Figure 
3  where  a  vehicle  and  a  tank  can  be  seen  closely  in  both  the  interpolated/sharpened 
version  and  the  super-resolved  version.  One  can  easily  see  the  clear  recovery  of  the 
edges  of  the  vehicle  and  tank  in  the  super-resolved  image  compared  to  the  artifactual 
versions  obtained  by  interpolation  and  sharpenning. 


(a)  (b) 


Figure  6:  Zoomed  areas:  (a)  Interpolated/sharpened,  (b)  super-resolved 


4  Conclusion 


A  method  of  multi-channel  super-resolution  has  been  proposed,  based  on  project¬ 
ing  image  frames  in  a  sequence  onto  a  more  optimal  Riesz  basis  and  exploiting  the 
temporal  bandwidth  to  increase  the  spatial  resolution.  Several  contributions  have 
been  made: 

•  A  unified  mathematical  framework  based  on  non-uniform  sampling  theory. 

•  Direct  extraction  of  channel  transfer  functions  from  the  input  data,  allowing  for 
the  design  and  implementation  of  sensor-dependent  reconstruction  and  adaptive 
estimation  of  variations  in  the  imaging  environment,  e.g.  the  PSF  and  the  noise. 

•  A  closed-form  solution  leading  to  a  high  speed  algorithm:  typically  15  seconds 
per  256x256  frame,  using  10-channel  data  fusion,  on  an  Ultra  Sparc  1. 

Thus  the  multi-frame  setup  commonly  used  in  the  literature  has  been  extended  to 
a  multi-channel  setup,  where  the  imaging  parameters  are  assumed  to  vary  from  one 
frame  to  another.  This  generalization  also  allows  extension  of  the  method  to  other 
sensor  types  such  as  radar  or  infra-red. 
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Appendix 


Proof  of  Theorem  2.2:  First,  a  well-known  result  is  applied  by  change  of  variables 
to  polar  coordinates: 

u  —  =  u>cos{(p)  V  —  Vk  =  cosin{(p) 

X  =  tcos(0)  y  =  tsin{0)  (21) 

This  leads  to  the  following  representation  of  the  Fourier  inversion  theorem  for  h: 

hk{t)exp{-itpk)  =  /  u>hk{oj)Jo{tu!)du>  (22) 

Jo 

where  pk  =  UkCos{6)  -b  Vksin{0),  hk{oj)  =  hk{u,v)  stands  for  the  Hankel  transform  of 
hk{t)  =  hk{x,y),  and  Jq  is  the  zero-order  Bessel  function  of  the  first  kind. 

Applying  the  Hankel  inversion  theorem  we  then  obtain 

hk{i*})  =  f  thk{t)exp{-itpk)Jo{tu;)dt  (23) 

Jo 

where  bk  is  the  radius  of  local  isotropy. 


Since  at  (u,u)  =  {uk,Vk)  we  have  hk{t>j)  =  0,  we  can  write 


f^k 

/  thk{t)exp{itpk)dt  =  0 

Jo 

(24) 

Substituting  from  (22)  into  this  last  integral  equation,  we  get 

r^k  roo  ^ 

/  t  u;hk{u;)Jo{t(jo)du;dt  =  0 

«/  0  V  0 

(25) 

which,  after  changing  the  order  of  integration,  yields 

yoo  _  rbk 

J  u)hk{<x>)  j  tJo{tu)dtduj  =  0 

(26) 

Using  the  well-known  identity 

f  tJo(t)dt  =  tJAat) 

Jo 

(27) 

we  deduce  that 

Jo  LO 

(28) 

Therefore 

/*co  ^ 

/  bkhk{u>)Ji{bkij^>)duj  =  0 
^0 

(29) 

Finally,  using  the  following  integral  relation  for  Bessel  functions  of  the  first  kind  [4]: 

Jn{buj)exp{±iu}a)  =  (a  -  -  fe2j  (30) 

we  conclude  that 

Yimhkhkipj)  =  exp{-iujtk)  (31) 

0/c  )-0 

where  tk  =  {x\  +  yl)2  and  {xk,yk)  is  a  constant  vector  in  T. 
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Errata 


1.  Page  7  paragraph  5  immediately  after  equation  (15)  should  read: 


where  hgi  =  and  hg2  =  /12  +  with  *  denoting  the  complex  conju¬ 

gate  and  I/ll  =  I/2I  =  |/|. 


2.  Page  7  equation  (16)  should  read: 


I  ^1^2  1=  c  I  I  where  c  =  |/| 
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